pub fn find_max_prime_factor(number: u128) -> u128 {
    if number < 2 {
        return number;
    }

    if is_prime(number) {
        return number;
    }

    let mut max_factor = 1;
    let mut stack = vec![number];

    while let Some(n) = stack.pop() {
        if n <= 1 {
            continue;
        }

        if is_prime(n) {
            if n > max_factor {
                max_factor = n;
            }
            continue;
        }

        let factor = pollards_rho(n);

        stack.push(factor);
        stack.push(n / factor);
    }

    max_factor
}

/// Miller-Rabin 素性测试
fn is_prime(n: u128) -> bool {
    if n < 2 {
        return false;
    }
    if n == 2 || n == 3 {
        return true;
    }
    if n % 2 == 0 {
        return false;
    }

    let d = n - 1;
    let s = d.trailing_zeros();
    let d = d >> s;

    // 选取一组确定性的测试 bases，适用于 u64
    let bases = [2u128, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];

    for &a in bases.iter() {
        if a >= n {
            break;
        }
        let mut x = modpow(a, d, n);
        if x == 1 || x == n - 1 {
            continue;
        }
        let mut continue_outer = false;
        for _ in 0..s - 1 {
            x = modmul(x, x, n);
            if x == n - 1 {
                continue_outer = true;
                break;
            }
        }
        if continue_outer {
            continue;
        }
        return false;
    }

    true
}

/// Pollard's Rho 算法
fn pollards_rho(n: u128) -> u128 {
    if n % 2 == 0 {
        return 2;
    }

    let mut x = 2u128;
    let mut y = 2u128;
    let mut c = 1u128;
    let mut d = 1u128;

    while d == 1 {
        x = (modmul(x, x, n) + c) % n;
        y = (modmul(y, y, n) + c) % n;
        y = (modmul(y, y, n) + c) % n;
        let diff = if x > y { x - y } else { y - x };
        d = gcd(diff, n);
        if d == n {
            c += 1;
            x = 2;
            y = 2;
            d = 1;
        }
    }

    d
}

fn modmul(a: u128, b: u128, m: u128) -> u128 {
    let mut result = 0u128;
    let mut a = a % m;
    let mut b = b;

    while b > 0 {
        if b % 2 == 1 {
            result = (result + a) % m;
        }
        a = (a * 2) % m;
        b /= 2;
    }

    result
}

fn modpow(mut base: u128, mut exponent: u128, modulus: u128) -> u128 {
    let mut result = 1u128;
    base = base % modulus;

    while exponent > 0 {
        if exponent % 2 == 1 {
            result = modmul(result, base, modulus);
        }
        base = modmul(base, base, modulus);
        exponent /= 2;
    }

    result
}

/// 计算最大公约数
fn gcd(mut a: u128, mut b: u128) -> u128 {
    while b != 0 {
        let tmp = a % b;
        a = b;
        b = tmp;
    }
    a
}